/*

  Simple cross-product, outputs in correct triangular form.

  - Coalescing memory access in the Column reads
  - Memory coalescing in the Row reads
  - 50% coalescing in the write (float8s)
  - All threads read in all data
  - Each thread works on a 2x2 square of data

 */

CUBE_KERNEL(shared2x2float4, cuFloatComplex *array_d, cuFloatComplex *product_d, 
	    const int write)
{
  CUBE_START;

  //get local thread ID
  int ty = threadIdx.y;
  int tx = threadIdx.x;
  int tid = ty*TILE_WIDTH + tx;

  //set frequency number from blockIdx.y
  int f = blockIdx.y;

  int Row, Col, blockX, blockY;
  CUBE_DEVICE_CALL(findPosition, Col, Row, blockX, blockY);
  Col*=2; Row*=2;

  //declare shared memory for input coalescing
  __shared__ float input[16*TILE_WIDTH]; // 4* for float4, 2* for 2x2 tile size

  //instantiate sum variables
  float sum11XXreal = 0.0, sum11XXimag = 0.0;
  float sum11XYreal = 0.0, sum11XYimag = 0.0;
  float sum11YXreal = 0.0, sum11YXimag = 0.0;
  float sum11YYreal = 0.0, sum11YYimag = 0.0;
  float sum12XXreal = 0.0, sum12XXimag = 0.0;
  float sum12XYreal = 0.0, sum12XYimag = 0.0;
  float sum12YXreal = 0.0, sum12YXimag = 0.0;
  float sum12YYreal = 0.0, sum12YYimag = 0.0;
  float sum21XXreal = 0.0, sum21XXimag = 0.0;
  float sum21XYreal = 0.0, sum21XYimag = 0.0;
  float sum21YXreal = 0.0, sum21YXimag = 0.0;
  float sum21YYreal = 0.0, sum21YYimag = 0.0;
  float sum22XXreal = 0.0, sum22XXimag = 0.0;
  float sum22XYreal = 0.0, sum22XYimag = 0.0;
  float sum22YXreal = 0.0, sum22YXimag = 0.0;
  float sum22YYreal = 0.0, sum22YYimag = 0.0;

#pragma unroll 4
  for(int t=0; t<NTIME; t++){
    __syncthreads();

    if(tid < 2*TILE_WIDTH){
      float4 temp = ((float4*)array_d)[t*NSTATION*NFREQUENCY + f*NSTATION + 2*blockX*TILE_WIDTH + tid];
      input[tid] = temp.x;
      input[tid + 2*TILE_WIDTH] = temp.y;
      input[tid + 4*TILE_WIDTH] = temp.z;
      input[tid + 6*TILE_WIDTH] = temp.w;
      CUBE_ADD_BYTES(16);
    } else if(tid >= warpSize && tid < warpSize+2*TILE_HEIGHT){
      float4 temp = ((float4*)array_d)[t*NSTATION*NFREQUENCY + f*NSTATION + 2*blockY*TILE_HEIGHT - warpSize + tid];
      input[tid - warpSize + 8*TILE_WIDTH] = temp.x;
      input[tid - warpSize + 2*TILE_WIDTH + 8*TILE_WIDTH] = temp.y;
      input[tid - warpSize + 4*TILE_WIDTH + 8*TILE_WIDTH] = temp.z;
      input[tid - warpSize + 6*TILE_WIDTH + 8*TILE_WIDTH] = temp.w;
      CUBE_ADD_BYTES(16);
    }

    __syncthreads();

    // read in shared data as individual floats to avoid bank conflicts
    float col1Xreal = input[2*tx];
    float col1Ximag = input[2*tx + 2*TILE_WIDTH];
    float col1Yreal = input[2*tx + 4*TILE_WIDTH];
    float col1Yimag = input[2*tx + 6*TILE_WIDTH];
    float col2Xreal = input[2*tx + 1];
    float col2Ximag = input[2*tx + 1 + 2*TILE_WIDTH];
    float col2Yreal = input[2*tx + 1 + 4*TILE_WIDTH];
    float col2Yimag = input[2*tx + 1 + 6*TILE_WIDTH];

    float row1Xreal = input[2*ty + 8*TILE_WIDTH];
    float row1Ximag = input[2*ty + 2*TILE_HEIGHT + 8*TILE_WIDTH];
    float row1Yreal = input[2*ty + 4*TILE_HEIGHT + 8*TILE_WIDTH];
    float row1Yimag = input[2*ty + 6*TILE_HEIGHT + 8*TILE_WIDTH];
    float row2Xreal = input[2*ty + 1 + 8*TILE_WIDTH];
    float row2Ximag = input[2*ty + 1 + 2*TILE_HEIGHT + 8*TILE_WIDTH];
    float row2Yreal = input[2*ty + 1 + 4*TILE_HEIGHT + 8*TILE_WIDTH];
    float row2Yimag = input[2*ty + 1 + 6*TILE_HEIGHT + 8*TILE_WIDTH];

    sum11XXreal += row1Xreal * col1Xreal;
    sum11XXreal += row1Ximag * col1Ximag;
    sum11XXimag += row1Ximag * col1Xreal;
    sum11XXimag -= row1Xreal * col1Ximag;
    sum11XYreal += row1Xreal * col1Yreal;
    sum11XYreal += row1Ximag * col1Yimag;
    sum11XYimag += row1Ximag * col1Yreal;
    sum11XYimag -= row1Xreal * col1Yimag;
    sum11YXreal += row1Yreal * col1Xreal;
    sum11YXreal += row1Yimag * col1Ximag;
    sum11YXimag += row1Yimag * col1Xreal;
    sum11YXimag -= row1Yreal * col1Ximag;
    sum11YYreal += row1Yreal * col1Yreal;
    sum11YYreal += row1Yimag * col1Yimag;
    sum11YYimag += row1Yimag * col1Yreal;
    sum11YYimag -= row1Yreal * col1Yimag;
    sum12XXreal += row1Xreal * col2Xreal;
    sum12XXreal += row1Ximag * col2Ximag;
    sum12XXimag += row1Ximag * col2Xreal;
    sum12XXimag -= row1Xreal * col2Ximag;
    sum12XYreal += row1Xreal * col2Yreal;
    sum12XYreal += row1Ximag * col2Yimag;
    sum12XYimag += row1Ximag * col2Yreal;
    sum12XYimag -= row1Xreal * col2Yimag;
    sum12YXreal += row1Yreal * col2Xreal;
    sum12YXreal += row1Yimag * col2Ximag;
    sum12YXimag += row1Yimag * col2Xreal;
    sum12YXimag -= row1Yreal * col2Ximag;
    sum12YYreal += row1Yreal * col2Yreal;
    sum12YYreal += row1Yimag * col2Yimag;
    sum12YYimag += row1Yimag * col2Yreal;
    sum12YYimag -= row1Yreal * col2Yimag;
    sum21XXreal += row2Xreal * col1Xreal;
    sum21XXreal += row2Ximag * col1Ximag;
    sum21XXimag += row2Ximag * col1Xreal;
    sum21XXimag -= row2Xreal * col1Ximag;
    sum21XYreal += row2Xreal * col1Yreal;
    sum21XYreal += row2Ximag * col1Yimag;
    sum21XYimag += row2Ximag * col1Yreal;
    sum21XYimag -= row2Xreal * col1Yimag;
    sum21YXreal += row2Yreal * col1Xreal;
    sum21YXreal += row2Yimag * col1Ximag;
    sum21YXimag += row2Yimag * col1Xreal;
    sum21YXimag -= row2Yreal * col1Ximag;
    sum21YYreal += row2Yreal * col1Yreal;
    sum21YYreal += row2Yimag * col1Yimag;
    sum21YYimag += row2Yimag * col1Yreal;
    sum21YYimag -= row2Yreal * col1Yimag;
    sum22XXreal += row2Xreal * col2Xreal;
    sum22XXreal += row2Ximag * col2Ximag;
    sum22XXimag += row2Ximag * col2Xreal;
    sum22XXimag -= row2Xreal * col2Ximag;
    sum22XYreal += row2Xreal * col2Yreal;
    sum22XYreal += row2Ximag * col2Yimag;
    sum22XYimag += row2Ximag * col2Yreal;
    sum22XYimag -= row2Xreal * col2Yimag;
    sum22YXreal += row2Yreal * col2Xreal;
    sum22YXreal += row2Yimag * col2Ximag;
    sum22YXimag += row2Yimag * col2Xreal;
    sum22YXimag -= row2Yreal * col2Ximag;
    sum22YYreal += row2Yreal * col2Yreal;
    sum22YYreal += row2Yimag * col2Yimag;
    sum22YYimag += row2Yimag * col2Yreal;
    sum22YYimag -= row2Yreal * col2Yimag;

    /*    sum11XXreal += row1Xreal * col1Xreal + row1Ximag * col1Ximag;
    sum11XXimag += row1Ximag * col1Xreal - row1Xreal * col1Ximag;
    sum11XYreal += row1Xreal * col1Yreal + row1Ximag * col1Yimag;
    sum11XYimag += row1Ximag * col1Yreal - row1Xreal * col1Yimag;
    sum11YXreal += row1Yreal * col1Xreal + row1Yimag * col1Ximag;
    sum11YXimag += row1Yimag * col1Xreal - row1Yreal * col1Ximag;
    sum11YYreal += row1Yreal * col1Yreal + row1Yimag * col1Yimag;
    sum11YYimag += row1Yimag * col1Yreal - row1Yreal * col1Yimag;

    sum12XXreal += row1Xreal * col2Xreal + row1Ximag * col2Ximag;
    sum12XXimag += row1Ximag * col2Xreal - row1Xreal * col2Ximag;
    sum12XYreal += row1Xreal * col2Yreal + row1Ximag * col2Yimag;
    sum12XYimag += row1Ximag * col2Yreal - row1Xreal * col2Yimag;
    sum12YXreal += row1Yreal * col2Xreal + row1Yimag * col2Ximag;
    sum12YXimag += row1Yimag * col2Xreal - row1Yreal * col2Ximag;
    sum12YYreal += row1Yreal * col2Yreal + row1Yimag * col2Yimag;
    sum12YYimag += row1Yimag * col2Yreal - row1Yreal * col2Yimag;

    sum21XXreal += row2Xreal * col1Xreal + row2Ximag * col1Ximag;
    sum21XXimag += row2Ximag * col1Xreal - row2Xreal * col1Ximag;
    sum21XYreal += row2Xreal * col1Yreal + row2Ximag * col1Yimag;
    sum21XYimag += row2Ximag * col1Yreal - row2Xreal * col1Yimag;
    sum21YXreal += row2Yreal * col1Xreal + row2Yimag * col1Ximag;
    sum21YXimag += row2Yimag * col1Xreal - row2Yreal * col1Ximag;
    sum21YYreal += row2Yreal * col1Yreal + row2Yimag * col1Yimag;
    sum21YYimag += row2Yimag * col1Yreal - row2Yreal * col1Yimag;

    sum22XXreal += row2Xreal * col2Xreal + row2Ximag * col2Ximag;
    sum22XXimag += row2Ximag * col2Xreal - row2Xreal * col2Ximag;
    sum22XYreal += row2Xreal * col2Yreal + row2Ximag * col2Yimag;
    sum22XYimag += row2Ximag * col2Yreal - row2Xreal * col2Yimag;
    sum22YXreal += row2Yreal * col2Xreal + row2Yimag * col2Ximag;
    sum22YXimag += row2Yimag * col2Xreal - row2Yreal * col2Ximag;
    sum22YYreal += row2Yreal * col2Yreal + row2Yimag * col2Yimag;
    sum22YYimag += row2Yimag * col2Yreal - row2Yreal * col2Yimag;*/
  } 


  if (write)
    CUBE_DEVICE_CALL(write2x2, Col, Row, product_d, sum11XXreal, sum11XXimag, sum11XYreal, 
		     sum11XYimag, sum11YXreal, sum11YXimag, sum11YYreal, sum11YYimag, 
		     sum12XXreal, sum12XXimag, sum12XYreal, sum12XYimag, sum12YXreal, 
		     sum12YXimag, sum12YYreal, sum12YYimag, sum21XXreal, sum21XXimag, 
		     sum21XYreal, sum21XYimag, sum21YXreal, sum21YXimag, sum21YYreal, 
		     sum21YYimag, sum22XXreal, sum22XXimag, sum22XYreal, sum22XYimag, 
		     sum22YXreal, sum22YXimag, sum22YYreal, sum22YYimag);

  if (Col < Row) {
    CUBE_ADD_BYTES(128);
    CUBE_ADD_FLOPS(NTIME*128);
  } else if (Col == Row) {
    CUBE_ADD_BYTES(96);
    CUBE_ADD_FLOPS(NTIME*96);
  }

  CUBE_END;
}
